pacman::p_load(readxl, httr, jsonlite, maptools, sf, sfdep, raster, spatstat, spNetwork, rgdal, sp, tmap, tidyverse)Project: 1st Order Spatial Point Patterns Analysis Methods
Geospatial Data Wrangling
For 1st order SPPA: airbnb_sf, mpsz_sf, sg_sf For 2nd order SPPA: same as above For NetSPPA: roadnetwork_sf For K function, K cross function, LCLQ: sevenele_sf, busstop_sf, hotel_sf, mrt_sf, mall_sf, attr_sf, unis_sf
read the 6 shp, kml, geojson files
busstop_sf <- st_read(dsn = "data/geospatial/busstop-022023-shp", layer="BusStop")Reading layer `BusStop' from data source
`C:\valtyl\IS415-GAA\Project\data\geospatial\busstop-022023-shp'
using driver `ESRI Shapefile'
Simple feature collection with 5159 features and 3 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 3970.122 ymin: 26482.1 xmax: 48284.56 ymax: 52983.82
Projected CRS: SVY21mpsz_sf <- st_read(dsn = "data/geospatial/MPSZ-2019-shp", layer="MPSZ-2019")Reading layer `MPSZ-2019' from data source
`C:\valtyl\IS415-GAA\Project\data\geospatial\MPSZ-2019-shp'
using driver `ESRI Shapefile'
Simple feature collection with 332 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 103.6057 ymin: 1.158699 xmax: 104.0885 ymax: 1.470775
Geodetic CRS: WGS 84network_sf <- st_read(dsn = "data/geospatial/roadnetwork-022023-shp", layer="RoadSectionLine")Reading layer `RoadSectionLine' from data source
`C:\valtyl\IS415-GAA\Project\data\geospatial\roadnetwork-022023-shp'
using driver `ESRI Shapefile'
Simple feature collection with 15300 features and 2 fields
Geometry type: MULTILINESTRING
Dimension: XY
Bounding box: xmin: 3248.394 ymin: 22755.1 xmax: 50043.6 ymax: 50135.95
Projected CRS: SVY21attr_sf <- st_read(dsn = "data/geospatial/tourism-shp", layer="TOURISM")Reading layer `TOURISM' from data source
`C:\valtyl\IS415-GAA\Project\data\geospatial\tourism-shp' using driver `ESRI Shapefile'
Simple feature collection with 109 features and 16 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 11380.23 ymin: 22858.67 xmax: 43659.54 ymax: 47627.69
Projected CRS: SVY21mrt_sf <- st_read(dsn = "data/geospatial/mrt-112022-shp", layer="Train_Station_Exit_Layer") Reading layer `Train_Station_Exit_Layer' from data source
`C:\valtyl\IS415-GAA\Project\data\geospatial\mrt-112022-shp'
using driver `ESRI Shapefile'
Simple feature collection with 562 features and 2 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 6134.086 ymin: 27499.7 xmax: 45356.36 ymax: 47865.92
Projected CRS: SVY21hotel_sf <- st_read("data/geospatial/hotels-2021-kml/hotel-locations.kml")Reading layer `HOTELS' from data source
`C:\valtyl\IS415-GAA\Project\data\geospatial\hotels-2021-kml\hotel-locations.kml'
using driver `KML'
Simple feature collection with 422 features and 2 fields
Geometry type: POINT
Dimension: XYZ
Bounding box: xmin: 103.6351 ymin: 1.245797 xmax: 103.9891 ymax: 1.419278
z_range: zmin: 0 zmax: 0
Geodetic CRS: WGS 84sg_sf <- st_read(dsn = "data/geospatial/costal-shp", layer="CostalOutline")Reading layer `CostalOutline' from data source
`C:\valtyl\IS415-GAA\Project\data\geospatial\costal-shp' using driver `ESRI Shapefile'
Simple feature collection with 60 features and 4 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: 2663.926 ymin: 16357.98 xmax: 56047.79 ymax: 50244.03
Projected CRS: SVY21read mall csv and convert to sf
mall_csv <- read_csv("data/geospatial/shoppingmalls-2019-csv/mall_coordinates_updated.csv")
mall_sf <- st_as_sf(mall_csv, coords = c("longitude", "latitude"), crs=4326)mall_sfSimple feature collection with 184 features and 2 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 103.6784 ymin: 1.263797 xmax: 103.9897 ymax: 1.448227
Geodetic CRS: WGS 84
# A tibble: 184 × 3
...1 name geometry
* <dbl> <chr> <POINT [°]>
1 0 100 AM (103.8435 1.274588)
2 1 112 KATONG (103.9051 1.305087)
3 2 313@SOMERSET (103.8377 1.301385)
4 3 321 CLEMENTI (103.765 1.312025)
5 4 600 @ TOA PAYOH (103.851 1.334042)
6 5 888 PLAZA (103.7953 1.437131)
7 6 AMK HUB (103.8485 1.369223)
8 7 ADMIRALTY PLACE (103.8018 1.439881)
9 8 ALEXANDRA RETAIL CENTRE (103.8014 1.273843)
10 9 ANCHORPOINT (103.8056 1.288935)
# … with 174 more rowsread airbnb csv and convert to sf
airbnb_csv <- read_csv("data/geospatial/airbnb-2022-csv/listings.csv")
airbnb_sf <- st_as_sf(airbnb_csv, coords = c("longitude", "latitude"), crs=4326)get_coords <- function(add_list){
# Create a data frame to store all retrieved coordinates
postal_coords <- data.frame()
# loop to go through each address in the list
for (i in add_list){
#print(i)
# response from OneMap API
r <- GET('https://developers.onemap.sg/commonapi/search?',
query=list(searchVal=i,
returnGeom='Y',
getAddrDetails='Y'))
data <- fromJSON(rawToChar(r$content))
found <- data$found
res <- data$results
# Create a new data frame for each address
new_row <- data.frame()
# If single result, append
if (found == 1){
postal <- res$POSTAL
lat <- res$LATITUDE
lng <- res$LONGITUDE
new_row <- data.frame(address= i, postal = postal, latitude = lat, longitude = lng)
}
# If multiple results, drop NIL and append top 1
else if (found > 1){
# Remove those with NIL as postal
res_sub <- res[res$POSTAL != "NIL", ]
# Set as NA first if no Postal
if (nrow(res_sub) == 0) {
new_row <- data.frame(address= i, postal = NA, latitude = NA, longitude = NA)
}
else{
top1 <- head(res_sub, n = 1)
postal <- top1$POSTAL
lat <- top1$LATITUDE
lng <- top1$LONGITUDE
new_row <- data.frame(address= i, postal = postal, latitude = lat, longitude = lng)
}
}
else {
new_row <- data.frame(address= i, postal = NA, latitude = NA, longitude = NA)
}
# Add the row to the main dataframe
postal_coords <- rbind(postal_coords, new_row)
}
return(postal_coords)
}read 7-11 xlsx and get coords
Show the code!
sevenele_xlsx <- read_xlsx("data/geospatial/711stores-032023-xlsx/7-11_convenience_stores.xlsx")
sevenele_list <- sort(unique(sevenele_xlsx$Postal))
sevenele_coords <- get_coords(sevenele_list)
sevenele_coords[(is.na(sevenele_coords$postal) | is.na(sevenele_coords$latitude) | is.na(sevenele_coords$longitude) | sevenele_coords$postal=="NIL"), ]fix 711 coords
Show the code!
sevenele_coords[64,]$postal <- "098975"
sevenele_coords[64,]$latitude <- "1.255169"
sevenele_coords[64,]$longitude <- "103.818508"
sevenele_coords[455,]$postal <- "828694"
sevenele_coords[455,]$latitude <- "1.420668"
sevenele_coords[455,]$longitude <- "103.912257"check if fixed
Show the code!
sevenele_coords[(is.na(sevenele_coords$postal) | is.na(sevenele_coords$latitude) | is.na(sevenele_coords$longitude) | sevenele_coords$postal=="NIL"), ]get 711 rds
Show the code!
sevenele_df <- left_join(sevenele_xlsx, sevenele_coords, by=c('Postal' = 'address'))
sevenele_rds <- write_rds(sevenele_df, "data/geospatial/711stores-032023-xlsx/sevenele_rds.rds")get 711 sf
seveneles <- read_rds("data/geospatial/711stores-032023-xlsx/sevenele_rds.rds")
sevenele_sf <- st_as_sf(seveneles,
coords = c("longitude",
"latitude"),
crs=4326) %>%
st_transform(crs = 3414)read universities and colleges xlsx and get coords
Show the code!
unis_xlsx <- read_xlsx("data/geospatial/universitiescolleges-2023-xlsx/universities_and_colleges.xlsx")
unis_list <- sort(unique(unis_xlsx$name))
unis_coords <- get_coords(unis_list)
unis_coords[(is.na(unis_coords$postal) | is.na(unis_coords$latitude) | is.na(unis_coords$longitude) | unis_coords$postal=="NIL"), ]fix uni coords
Show the code!
unis_coords[1,]$postal <- "179104"
unis_coords[1,]$latitude <- "1.29006"
unis_coords[1,]$longitude <- "103.8497"
unis_coords[2,]$postal <- "229469"
unis_coords[2,]$latitude <- "1.31648392303"
unis_coords[2,]$longitude <- "103.873964955"
unis_coords[4,]$postal <- "408601"
unis_coords[4,]$latitude <- "1.3196562"
unis_coords[4,]$longitude <- "103.8923565"
unis_coords[5,]$postal <- "059405"
unis_coords[5,]$latitude <- "1.28799"
unis_coords[5,]$longitude <- "103.84692"
unis_coords[6,]$postal <- "188946"
unis_coords[6,]$latitude <- "1.3462717"
unis_coords[6,]$longitude <- "103.95302009"
unis_coords[7,]$postal <- "329162"
unis_coords[7,]$latitude <- "1.3270276324741"
unis_coords[7,]$longitude <- "103.85336279869"
unis_coords[10,]$postal <- "188647"
unis_coords[10,]$latitude <- "1.2990616"
unis_coords[10,]$longitude <- "103.8494881"
unis_coords[11,]$postal <- "088383"
unis_coords[11,]$latitude <- "1.28013"
unis_coords[11,]$longitude <- "103.84193"
unis_coords[12,]$postal <- "049319"
unis_coords[12,]$latitude <- "1.28469"
unis_coords[12,]$longitude <- "103.85269"
unis_coords[14,]$postal <- "138676"
unis_coords[14,]$latitude <- "1.299122"
unis_coords[14,]$longitude <- "103.788025"
unis_coords[16,]$postal <- "228095"
unis_coords[16,]$latitude <- "1.30209070954"
unis_coords[16,]$longitude <- "103.849560613"
unis_coords[17,]$postal <- "188306"
unis_coords[17,]$latitude <- "1.3001799"
unis_coords[17,]$longitude <- "103.84924"
unis_coords[19,]$postal <- "069542"
unis_coords[19,]$latitude <- "1.2750794"
unis_coords[19,]$longitude <- "103.8462627"
unis_coords[20,]$postal <- "148951"
unis_coords[20,]$latitude <- "1.2970726"
unis_coords[20,]$longitude <- "103.8011673"
unis_coords[22,]$postal <- "189655"
unis_coords[22,]$latitude <- "1.30020700501"
unis_coords[22,]$longitude <- "103.851710836"
unis_coords[26,]$postal <- "248922"
unis_coords[26,]$latitude <- "1.31532"
unis_coords[26,]$longitude <- "103.77994"
unis_coords[30,]$postal <- "599491"
unis_coords[30,]$latitude <- "1.326"
unis_coords[30,]$longitude <- "103.79999"
unis_coords[37,]$postal <- "550266"
unis_coords[37,]$latitude <- "1.35324"
unis_coords[37,]$longitude <- "103.87145"
unis_coords[38,]$postal <- "139660"
unis_coords[38,]$latitude <- "1.3079"
unis_coords[38,]$longitude <- "103.777"check fix
Show the code!
unis_coords[(is.na(unis_coords$postal) | is.na(unis_coords$latitude) | is.na(unis_coords$longitude) | unis_coords$postal=="NIL"), ]get uni rds
Show the code!
unis_df <- left_join(unis_xlsx, unis_coords, by=c('name' = 'address'))
unis_rds <- write_rds(unis_df, "data/geospatial/universitiescolleges-2023-xlsx/unis_rds.rds")get unis sf
unis <- read_rds("data/geospatial/universitiescolleges-2023-xlsx/unis_rds.rds")
unis_sf <- st_as_sf(unis,
coords = c("longitude",
"latitude"),
crs=4326) %>%
st_transform(crs = 3414)remove unnecessary columns (except mpsz_sf, network_sf)
attr_sf <- attr_sf %>% select(c(5))
busstop_sf <- busstop_sf %>% select(c(1))
hotel_sf <- hotel_sf %>% select(c(1))
mall_sf <- mall_sf %>% select(c(2))
mrt_sf <- mrt_sf %>% select(c(1))
sevenele_sf <- sevenele_sf %>% select(c(1))
unis_sf <- unis_sf %>% select(c(1))
airbnb_sf <- airbnb_sf %>% select(c(2,6,7,8))remove z dimensions
hotel_sf <- st_zm(hotel_sf)check for invalid geometries
length(which(st_is_valid(attr_sf) == FALSE))[1] 0length(which(st_is_valid(busstop_sf) == FALSE))[1] 0length(which(st_is_valid(hotel_sf) == FALSE))[1] 0length(which(st_is_valid(mall_sf) == FALSE))[1] 0length(which(st_is_valid(mrt_sf) == FALSE))[1] 0length(which(st_is_valid(sevenele_sf) == FALSE))[1] 0length(which(st_is_valid(unis_sf) == FALSE))[1] 0length(which(st_is_valid(mpsz_sf) == FALSE))[1] 6length(which(st_is_valid(network_sf) == FALSE))[1] 0length(which(st_is_valid(airbnb_sf) == FALSE))[1] 0length(which(st_is_valid(sg_sf) == FALSE))[1] 1hotel_sf <- st_zm(hotel_sf)fix mpsz invalid geometry and check again
mpsz_sf <- st_make_valid(mpsz_sf)
length(which(st_is_valid(mpsz_sf) == FALSE))[1] 0fix sg invalid geometry and check again
sg_sf <- st_make_valid(sg_sf)
length(which(st_is_valid(sg_sf) == FALSE))[1] 0dealing w missing values
attr_sf[rowSums(is.na(attr_sf))!=0,]Simple feature collection with 0 features and 1 field
Bounding box: xmin: NA ymin: NA xmax: NA ymax: NA
Projected CRS: SVY21
[1] PAGETITLE geometry
<0 rows> (or 0-length row.names)busstop_sf[rowSums(is.na(busstop_sf))!=0,]Simple feature collection with 0 features and 1 field
Bounding box: xmin: NA ymin: NA xmax: NA ymax: NA
Projected CRS: SVY21
[1] BUS_STOP_N geometry
<0 rows> (or 0-length row.names)hotel_sf[rowSums(is.na(hotel_sf))!=0,]Simple feature collection with 0 features and 1 field
Bounding box: xmin: NA ymin: NA xmax: NA ymax: NA
Geodetic CRS: WGS 84
[1] Name geometry
<0 rows> (or 0-length row.names)mall_sf[rowSums(is.na(mall_sf))!=0,]Simple feature collection with 0 features and 1 field
Bounding box: xmin: NA ymin: NA xmax: NA ymax: NA
Geodetic CRS: WGS 84
# A tibble: 0 × 2
# … with 2 variables: name <chr>, geometry <GEOMETRY [°]>mrt_sf[rowSums(is.na(mrt_sf))!=0,]Simple feature collection with 0 features and 1 field
Bounding box: xmin: NA ymin: NA xmax: NA ymax: NA
Projected CRS: SVY21
[1] stn_name geometry
<0 rows> (or 0-length row.names)sevenele_sf[rowSums(is.na(sevenele_sf))!=0,]Simple feature collection with 0 features and 1 field
Bounding box: xmin: NA ymin: NA xmax: NA ymax: NA
Projected CRS: SVY21 / Singapore TM
# A tibble: 0 × 2
# … with 2 variables: Name <chr>, geometry <GEOMETRY [m]>unis_sf[rowSums(is.na(unis_sf))!=0,]Simple feature collection with 0 features and 1 field
Bounding box: xmin: NA ymin: NA xmax: NA ymax: NA
Projected CRS: SVY21 / Singapore TM
# A tibble: 0 × 2
# … with 2 variables: name <chr>, geometry <GEOMETRY [m]>mpsz_sf[rowSums(is.na(mpsz_sf))!=0,]Simple feature collection with 0 features and 6 fields
Bounding box: xmin: NA ymin: NA xmax: NA ymax: NA
Geodetic CRS: WGS 84
[1] SUBZONE_N SUBZONE_C PLN_AREA_N PLN_AREA_C REGION_N REGION_C geometry
<0 rows> (or 0-length row.names)network_sf[rowSums(is.na(network_sf))!=0,]Simple feature collection with 15300 features and 2 fields
Geometry type: MULTILINESTRING
Dimension: XY
Bounding box: xmin: 3248.394 ymin: 22755.1 xmax: 50043.6 ymax: 50135.95
Projected CRS: SVY21
First 10 features:
RD_CD RD_CD_DESC geometry
1 <NA> SEA BREEZE ROAD MULTILINESTRING ((41260.55 ...
2 <NA> HILLVIEW AVENUE MULTILINESTRING ((20404.4 3...
3 <NA> SERANGOON CENTRAL MULTILINESTRING ((32508.6 3...
4 <NA> CHANGI SOUTH STREET 3 MULTILINESTRING ((42949.75 ...
5 <NA> PANDAN ROAD MULTILINESTRING ((18683.23 ...
6 <NA> WEST COAST GROVE MULTILINESTRING ((19396.72 ...
7 <NA> NEYTHAL ROAD MULTILINESTRING ((12917.95 ...
8 <NA> WILLOW AVENUE MULTILINESTRING ((32638.54 ...
9 <NA> CHOA CHU KANG NORTH 5 MULTILINESTRING ((18182.78 ...
10 <NA> CASSIA CRESCENT MULTILINESTRING ((33711.34 ...airbnb_sf[rowSums(is.na(airbnb_sf))!=0,]Simple feature collection with 0 features and 4 fields
Bounding box: xmin: NA ymin: NA xmax: NA ymax: NA
Geodetic CRS: WGS 84
# A tibble: 0 × 5
# … with 5 variables: name <chr>, neighbourhood <chr>, room_type <chr>,
# price <dbl>, geometry <GEOMETRY [°]>sg_sf[rowSums(is.na(sg_sf))!=0,]Simple feature collection with 0 features and 4 fields
Bounding box: xmin: NA ymin: NA xmax: NA ymax: NA
Projected CRS: SVY21
[1] GDO_GID MSLINK MAPID COSTAL_NAM geometry
<0 rows> (or 0-length row.names)verify coordinate system
st_crs(attr_sf)Coordinate Reference System:
User input: SVY21
wkt:
PROJCRS["SVY21",
BASEGEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ID["EPSG",6326]],
PRIMEM["Greenwich",0,
ANGLEUNIT["Degree",0.0174532925199433]]],
CONVERSION["unnamed",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",1.36666666666667,
ANGLEUNIT["Degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",103.833333333333,
ANGLEUNIT["Degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",1,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",28001.642,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",38744.572,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]]]st_crs(busstop_sf)Coordinate Reference System:
User input: SVY21
wkt:
PROJCRS["SVY21",
BASEGEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ID["EPSG",6326]],
PRIMEM["Greenwich",0,
ANGLEUNIT["Degree",0.0174532925199433]]],
CONVERSION["unnamed",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",1.36666666666667,
ANGLEUNIT["Degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",103.833333333333,
ANGLEUNIT["Degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",1,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",28001.642,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",38744.572,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]]]st_crs(hotel_sf)Coordinate Reference System:
User input: WGS 84
wkt:
GEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]]st_crs(mall_sf)Coordinate Reference System:
User input: EPSG:4326
wkt:
GEOGCRS["WGS 84",
ENSEMBLE["World Geodetic System 1984 ensemble",
MEMBER["World Geodetic System 1984 (Transit)"],
MEMBER["World Geodetic System 1984 (G730)"],
MEMBER["World Geodetic System 1984 (G873)"],
MEMBER["World Geodetic System 1984 (G1150)"],
MEMBER["World Geodetic System 1984 (G1674)"],
MEMBER["World Geodetic System 1984 (G1762)"],
MEMBER["World Geodetic System 1984 (G2139)"],
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[2.0]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
USAGE[
SCOPE["Horizontal component of 3D system."],
AREA["World."],
BBOX[-90,-180,90,180]],
ID["EPSG",4326]]st_crs(mrt_sf)Coordinate Reference System:
User input: SVY21
wkt:
PROJCRS["SVY21",
BASEGEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ID["EPSG",6326]],
PRIMEM["Greenwich",0,
ANGLEUNIT["Degree",0.0174532925199433]]],
CONVERSION["unnamed",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",1.36666666666667,
ANGLEUNIT["Degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",103.833333333333,
ANGLEUNIT["Degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",1,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",28001.642,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",38744.572,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]]]st_crs(sevenele_sf)Coordinate Reference System:
User input: EPSG:3414
wkt:
PROJCRS["SVY21 / Singapore TM",
BASEGEOGCRS["SVY21",
DATUM["SVY21",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4757]],
CONVERSION["Singapore Transverse Mercator",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",1.36666666666667,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",103.833333333333,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",1,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",28001.642,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",38744.572,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["northing (N)",north,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["easting (E)",east,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Cadastre, engineering survey, topographic mapping."],
AREA["Singapore - onshore and offshore."],
BBOX[1.13,103.59,1.47,104.07]],
ID["EPSG",3414]]st_crs(unis_sf)Coordinate Reference System:
User input: EPSG:3414
wkt:
PROJCRS["SVY21 / Singapore TM",
BASEGEOGCRS["SVY21",
DATUM["SVY21",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4757]],
CONVERSION["Singapore Transverse Mercator",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",1.36666666666667,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",103.833333333333,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",1,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",28001.642,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",38744.572,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["northing (N)",north,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["easting (E)",east,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Cadastre, engineering survey, topographic mapping."],
AREA["Singapore - onshore and offshore."],
BBOX[1.13,103.59,1.47,104.07]],
ID["EPSG",3414]]st_crs(mpsz_sf)Coordinate Reference System:
User input: WGS 84
wkt:
GEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["latitude",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["longitude",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]]st_crs(network_sf)Coordinate Reference System:
User input: SVY21
wkt:
PROJCRS["SVY21",
BASEGEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ID["EPSG",6326]],
PRIMEM["Greenwich",0,
ANGLEUNIT["Degree",0.0174532925199433]]],
CONVERSION["unnamed",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",1.36666666666667,
ANGLEUNIT["Degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",103.833333333333,
ANGLEUNIT["Degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",1,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",28001.642,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",38744.572,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]]]st_crs(airbnb_sf)Coordinate Reference System:
User input: EPSG:4326
wkt:
GEOGCRS["WGS 84",
ENSEMBLE["World Geodetic System 1984 ensemble",
MEMBER["World Geodetic System 1984 (Transit)"],
MEMBER["World Geodetic System 1984 (G730)"],
MEMBER["World Geodetic System 1984 (G873)"],
MEMBER["World Geodetic System 1984 (G1150)"],
MEMBER["World Geodetic System 1984 (G1674)"],
MEMBER["World Geodetic System 1984 (G1762)"],
MEMBER["World Geodetic System 1984 (G2139)"],
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[2.0]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
USAGE[
SCOPE["Horizontal component of 3D system."],
AREA["World."],
BBOX[-90,-180,90,180]],
ID["EPSG",4326]]st_crs(sg_sf)Coordinate Reference System:
User input: SVY21
wkt:
PROJCRS["SVY21",
BASEGEOGCRS["SVY21[WGS84]",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ID["EPSG",6326]],
PRIMEM["Greenwich",0,
ANGLEUNIT["Degree",0.0174532925199433]]],
CONVERSION["unnamed",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",1.36666666666667,
ANGLEUNIT["Degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",103.833333333333,
ANGLEUNIT["Degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",1,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",28001.642,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",38744.572,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]]]transform coordinate system
# with st_set_crs(), we can assign the appropriate ESPG Code
attr_sf <- st_set_crs(attr_sf, 3414)
busstop_sf <- st_set_crs(busstop_sf, 3414)
mrt_sf <- st_set_crs(mrt_sf, 3414)
network_sf <- st_set_crs(network_sf, 3414)
sg_sf <- st_set_crs(sg_sf, 3414)
# with st_transform(), we can change from one CRS to another
mpsz_sf <- st_transform(mpsz_sf, crs=3414)
hotel_sf <- st_transform(hotel_sf, crs=3414)
mall_sf <- st_transform(mall_sf, crs=3414)
airbnb_sf <- st_transform(airbnb_sf, crs=3414)
# sevenele_sf and unis_sf are in the correct CRS and ESPG codePreparing for 1st Order SPPA
convert sf to sp’s spatial class
airbnb <- as_Spatial(airbnb_sf)
mpsz <- as_Spatial(mpsz_sf)
sg <- as_Spatial(sg_sf)convert sp’s spatial into generic sp format
airbnb_sp <- as(airbnb, "SpatialPoints")
sg_sp <- as(sg, "SpatialPolygons")airbnbclass : SpatialPointsDataFrame
features : 3037
extent : 7406.989, 44064, 25352.92, 48321.84 (xmin, xmax, ymin, ymax)
crs : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
variables : 4
names : name, neighbourhood, room_type, price
min values : -Master bedroom with attached bathroom at CBD, Ang Mo Kio, Entire home/apt, 0
max values : 阿裕尼公寓普通房急需女搭房一名联系人:郑小姐 电话：, Yishun, Shared room, 135140 convert generic sp format into spatstat’s ppp format
airbnb_ppp <- as(airbnb_sp, "ppp")
airbnb_pppPlanar point pattern: 3037 points
window: rectangle = [7406.99, 44064] x [25352.92, 48321.84] unitscheck number of duplicate point events
sum(multiplicity(airbnb_ppp) > 1)[1] 211handle duplicate points:
airbnb_ppp_jit <- rjitter(airbnb_ppp,
retry=TRUE,
nsim=1,
drop=TRUE)check duplicated points:
any(duplicated(airbnb_ppp_jit))[1] FALSEcreate owin
# sg_owin <- as(sg_sp, "owin")
central = mpsz_sf[mpsz_sf$REGION_N=='CENTRAL REGION',]
central_owin <- central %>%
as ('Spatial') %>%
as ('SpatialPolygons') %>%
as ('owin')
plot(central_owin)
combine airbnb point events obj and owin obj
airbnbSG_ppp = airbnb_ppp_jit[central_owin]summary
summary(airbnbSG_ppp)Planar point pattern: 2476 points
Average intensity 1.815334e-05 points per square unit
Coordinates are given to 3 decimal places
i.e. rounded to the nearest multiple of 0.001 units
Window: polygonal boundary
164 separate polygons (14 holes)
vertices area relative.area
polygon 1 299 1.84404e+06 1.35e-02
polygon 2 164 3.92563e+05 2.88e-03
polygon 3 238 5.06591e+05 3.71e-03
polygon 4 132 3.88733e+05 2.85e-03
polygon 5 255 1.59034e+06 1.17e-02
polygon 6 1020 1.27781e+06 9.37e-03
polygon 7 (hole) 4 -1.23063e-04 -9.02e-13
polygon 8 (hole) 3 -1.06327e-01 -7.80e-10
polygon 9 (hole) 5 -5.18342e-01 -3.80e-09
polygon 10 (hole) 3 -2.53510e-02 -1.86e-10
polygon 11 (hole) 4 -7.56223e-02 -5.54e-10
polygon 12 (hole) 4 -7.19058e-01 -5.27e-09
polygon 13 211 4.70522e+05 3.45e-03
polygon 14 155 2.67503e+05 1.96e-03
polygon 15 129 9.53762e+04 6.99e-04
polygon 16 94 5.96175e+04 4.37e-04
polygon 17 59 3.43163e+04 2.52e-04
polygon 18 10 4.90997e+02 3.60e-06
polygon 19 6 4.50376e+02 3.30e-06
polygon 20 4 2.69408e+02 1.98e-06
polygon 21 1424 4.87153e+06 3.57e-02
polygon 22 (hole) 4 -2.46124e-04 -1.80e-12
polygon 23 (hole) 3 -3.69186e-03 -2.71e-11
polygon 24 (hole) 3 -2.73199e-02 -2.00e-10
polygon 25 (hole) 11 -8.36614e+01 -6.13e-07
polygon 26 75 1.73525e+04 1.27e-04
polygon 27 40 1.38603e+04 1.02e-04
polygon 28 83 5.28926e+03 3.88e-05
polygon 29 137 3.22285e+03 2.36e-05
polygon 30 147 3.10396e+03 2.28e-05
polygon 31 106 3.04131e+03 2.23e-05
polygon 32 45 2.51228e+03 1.84e-05
polygon 33 438 3.45157e+06 2.53e-02
polygon 34 83 1.03237e+05 7.57e-04
polygon 35 102 1.12730e+06 8.27e-03
polygon 36 829 2.70003e+06 1.98e-02
polygon 37 (hole) 23 -1.99689e+01 -1.46e-07
polygon 38 (hole) 35 -1.38430e+02 -1.01e-06
polygon 39 (hole) 19 -4.37731e+00 -3.21e-08
polygon 40 (hole) 270 -1.21366e+03 -8.90e-06
polygon 41 88 5.33462e+05 3.91e-03
polygon 42 95 1.45518e+05 1.07e-03
polygon 43 55 6.35706e+05 4.66e-03
polygon 44 53 2.76828e+05 2.03e-03
polygon 45 114 6.36655e+04 4.67e-04
polygon 46 83 1.96620e+05 1.44e-03
polygon 47 33 3.65334e+05 2.68e-03
polygon 48 110 1.45503e+06 1.07e-02
polygon 49 135 8.53202e+05 6.26e-03
polygon 50 196 1.07072e+06 7.85e-03
polygon 51 47 5.33012e+05 3.91e-03
polygon 52 85 4.42296e+05 3.24e-03
polygon 53 38 4.11723e+05 3.02e-03
polygon 54 227 5.87221e+05 4.31e-03
polygon 55 35 3.94370e+04 2.89e-04
polygon 56 96 1.88768e+05 1.38e-03
polygon 57 59 1.33007e+05 9.75e-04
polygon 58 46 4.48128e+05 3.29e-03
polygon 59 31 5.21199e+05 3.82e-03
polygon 60 17 3.50788e+05 2.57e-03
polygon 61 54 2.61841e+05 1.92e-03
polygon 62 152 1.63038e+06 1.20e-02
polygon 63 180 5.61608e+05 4.12e-03
polygon 64 47 1.60809e+05 1.18e-03
polygon 65 49 5.95241e+05 4.36e-03
polygon 66 49 3.87611e+05 2.84e-03
polygon 67 59 1.03038e+06 7.55e-03
polygon 68 83 5.51731e+05 4.05e-03
polygon 69 69 2.90185e+05 2.13e-03
polygon 70 217 1.08479e+06 7.95e-03
polygon 71 41 6.28892e+05 4.61e-03
polygon 72 226 1.82685e+06 1.34e-02
polygon 73 55 2.93699e+05 2.15e-03
polygon 74 247 5.57280e+05 4.09e-03
polygon 75 48 5.56817e+04 4.08e-04
polygon 76 60 1.16331e+05 8.53e-04
polygon 77 343 2.00345e+06 1.47e-02
polygon 78 129 2.43459e+06 1.78e-02
polygon 79 58 3.10535e+05 2.28e-03
polygon 80 114 1.38034e+06 1.01e-02
polygon 81 134 1.95176e+06 1.43e-02
polygon 82 270 4.51538e+05 3.31e-03
polygon 83 123 1.72655e+05 1.27e-03
polygon 84 137 3.36218e+05 2.47e-03
polygon 85 62 7.42203e+05 5.44e-03
polygon 86 114 1.07899e+06 7.91e-03
polygon 87 319 4.60551e+05 3.38e-03
polygon 88 198 5.43484e+05 3.98e-03
polygon 89 51 2.78303e+05 2.04e-03
polygon 90 297 8.86957e+05 6.50e-03
polygon 91 182 2.18809e+05 1.60e-03
polygon 92 129 2.00787e+05 1.47e-03
polygon 93 169 7.10568e+05 5.21e-03
polygon 94 34 7.48683e+05 5.49e-03
polygon 95 173 3.68483e+05 2.70e-03
polygon 96 294 7.60622e+06 5.58e-02
polygon 97 238 2.21999e+05 1.63e-03
polygon 98 130 2.80176e+05 2.05e-03
polygon 99 140 2.14251e+05 1.57e-03
polygon 100 83 1.73123e+05 1.27e-03
polygon 101 192 5.91779e+05 4.34e-03
polygon 102 174 1.75348e+06 1.29e-02
polygon 103 192 3.40742e+05 2.50e-03
polygon 104 218 3.29438e+05 2.42e-03
polygon 105 87 1.70664e+05 1.25e-03
polygon 106 217 1.34616e+06 9.87e-03
polygon 107 27 1.71336e+05 1.26e-03
polygon 108 84 4.96255e+04 3.64e-04
polygon 109 198 1.93991e+05 1.42e-03
polygon 110 77 1.20171e+05 8.81e-04
polygon 111 273 6.14923e+05 4.51e-03
polygon 112 108 1.02656e+06 7.53e-03
polygon 113 154 1.67538e+05 1.23e-03
polygon 114 81 1.16002e+06 8.50e-03
polygon 115 86 9.63497e+05 7.06e-03
polygon 116 103 1.42508e+06 1.04e-02
polygon 117 74 1.22990e+06 9.02e-03
polygon 118 75 1.26341e+06 9.26e-03
polygon 119 50 3.69771e+05 2.71e-03
polygon 120 96 1.10727e+06 8.12e-03
polygon 121 126 3.38725e+06 2.48e-02
polygon 122 94 1.87800e+06 1.38e-02
polygon 123 40 8.67749e+05 6.36e-03
polygon 124 55 6.39144e+05 4.69e-03
polygon 125 54 4.11402e+05 3.02e-03
polygon 126 75 4.18656e+05 3.07e-03
polygon 127 74 2.18107e+05 1.60e-03
polygon 128 104 2.13519e+05 1.57e-03
polygon 129 119 4.85574e+05 3.56e-03
polygon 130 102 7.56923e+05 5.55e-03
polygon 131 117 3.51306e+05 2.58e-03
polygon 132 67 9.47451e+05 6.95e-03
polygon 133 93 1.05203e+06 7.71e-03
polygon 134 92 4.10940e+05 3.01e-03
polygon 135 72 8.39679e+05 6.16e-03
polygon 136 184 1.22706e+06 9.00e-03
polygon 137 120 5.58761e+05 4.10e-03
polygon 138 212 2.09608e+06 1.54e-02
polygon 139 119 2.15829e+06 1.58e-02
polygon 140 140 1.34746e+06 9.88e-03
polygon 141 381 2.03851e+06 1.49e-02
polygon 142 125 2.57843e+06 1.89e-02
polygon 143 112 7.35502e+05 5.39e-03
polygon 144 132 1.31911e+06 9.67e-03
polygon 145 122 1.37683e+06 1.01e-02
polygon 146 129 1.92662e+06 1.41e-02
polygon 147 74 1.79446e+06 1.32e-02
polygon 148 136 3.14493e+06 2.31e-02
polygon 149 195 2.63648e+06 1.93e-02
polygon 150 80 1.05717e+06 7.75e-03
polygon 151 69 4.39648e+05 3.22e-03
polygon 152 61 4.46241e+05 3.27e-03
polygon 153 72 5.72500e+05 4.20e-03
polygon 154 79 8.13382e+05 5.96e-03
polygon 155 97 1.03728e+06 7.61e-03
polygon 156 182 2.77329e+06 2.03e-02
polygon 157 131 3.46851e+06 2.54e-02
polygon 158 71 2.35998e+05 1.73e-03
polygon 159 244 4.70674e+05 3.45e-03
polygon 160 97 9.13532e+04 6.70e-04
polygon 161 128 1.71195e+06 1.26e-02
polygon 162 163 2.96147e+06 2.17e-02
polygon 163 268 3.84951e+06 2.82e-02
polygon 164 102 1.96415e+06 1.44e-02
enclosing rectangle: [18752.4, 37592.19] x [21678.34, 38891.36] units
(18840 x 17210 units)
Window area = 136394000 square units
Fraction of frame area: 0.4211st Order SPPA
rescale to convert from m to km
airbnbSG_ppp.km <- rescale(airbnbSG_ppp, 1000, "km")Kernel Density Estimation (automatic bandwidth)
kde_airbnbSG.bw <- density(airbnbSG_ppp.km,
sigma=bw.diggle,
edge=TRUE,
kernel="gaussian") other methods to do: bw.CvL(), bw.scott(), bw.ppl() other smoothing kernels: epanechnikov, quartic, disc
plot(kde_airbnbSG.bw)
KDE (fixed bandwidth)
kde_airbnbSG_600 <- density(airbnbSG_ppp.km, sigma=0.6, edge=TRUE, kernel="gaussian") #0.6 sigma = 600 m for bandwidth. 600m = 0.6km
plot(kde_airbnbSG_600)
sigma: range from 0.1 to 10 different kernels
KDE (adaptive bandwidth)
kde_airbnbSG_adaptive <- adaptive.density(airbnbSG_ppp.km, method="kernel")
plot(kde_airbnbSG_adaptive)
convert KDE output into grid (do for all)
gridded_kde_airbnbSG_bw <- as.SpatialGridDataFrame.im(kde_airbnbSG.bw)
# convert grid output into raster
kde_airbnbSG_bw_raster <- raster(gridded_kde_airbnbSG_bw)
# assign projection systems
projection(kde_airbnbSG_bw_raster) <- CRS("+init=EPSG:3414 +datum=WGS84 +units=km")
kde_airbnbSG_bw_rasterclass : RasterLayer
dimensions : 128, 128, 16384 (nrow, ncol, ncell)
resolution : 0.1471859, 0.1344767 (x, y)
extent : 18.7524, 37.59219, 21.67834, 38.89136 (xmin, xmax, ymin, ymax)
crs : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +datum=WGS84 +units=km +no_defs
source : memory
names : v
values : -3.377976e-13, 2071.377 (min, max)kl = mpsz_sf[mpsz_sf$PLN_AREA_N == "KALLANG",]
dc = mpsz_sf[mpsz_sf$PLN_AREA_N == "DOWNTOWN CORE",]
ot = mpsz_sf[mpsz_sf$PLN_AREA_N == "OUTRAM",]
rc = mpsz_sf[mpsz_sf$PLN_AREA_N == "ROCHOR",]mapex <- st_bbox(central)final visualisation tmap
tmap_mode("plot")
tm_shape(kde_airbnbSG_bw_raster, bbox=mapex) +
tm_raster("v") +
tm_layout(legend.position = c("right", "bottom"), frame = FALSE)
2nd Order SPPA G/F/K/L Function
G Function estimation
G_airbnb = Gest(airbnbSG_ppp, correction =c("rs", "Hanisch"))
plot(G_airbnb, xlim=c(0,500))other correction: “best” Gest() - corrections: “none”, “rs”, “km”, “Hanisch”, “best”, “all” Fest() - corrections: “none”, “rs”, “km”, “cs”, “best”, “all” Kest() - corrections: “none”, “border”, “bord.modif”, “isotropic”, “Ripley”, “translate”, “translation”, “rigid”, “none”, “good”, “best”, “all” Lest() - same as Kest()
Complete Spatial Randomness Test
G_airbnb.csr <- envelope(airbnbSG_ppp, Gest, nsim = 999)
# G_airbnb.csr <- envelope(airbnbSG_ppp, Gest, correction="all", nsim = 999)nsim: number btwn 0 to 1000
plot to show
plot(G_airbnb.csr)K Function estimation
K_airbnb = Kest(airbnbSG_ppp, correction =c("best"))
plot(K_airbnb, xlim=c(0,500))Complete Spatial Randomness Test
K_airbnb.csr <- envelope(airbnbSG_ppp, Kest, nsim = 999)
# G_airbnb.csr <- envelope(airbnbSG_ppp, Gest, correction="all", nsim = 999)nsim: number btwn 0 to 1000
plot to show
plot(K_airbnb.csr)visualisation
tmap_mode("plot")
tm_shape(mpsz_sf) +
tm_borders(alpha = 0.5) +
tm_shape(airbnb_sf) +
tm_dots(col = "red", size = 0.05) +
tm_layout(main.title = "Airbnbs",
main.title.position = 'center',
main.title.size = 0.8,
frame = TRUE)
occurrence <- airbnb_sf %>% count(airbnb_sf$neighbourhood)occurrence <- occurrence[order(-occurrence$n),]occurrenceSimple feature collection with 43 features and 2 fields
Geometry type: GEOMETRY
Dimension: XY
Bounding box: xmin: 7406.989 ymin: 25352.92 xmax: 44064 ymax: 48321.84
Projected CRS: SVY21 / Singapore TM
# A tibble: 43 × 3
`airbnb_sf$neighbourhood` n geometry
<chr> <int> <MULTIPOINT [m]>
1 Kallang 353 ((29508.14 32976.25), (29791.93 33174.18), (…
2 Downtown Core 325 ((28945.03 28564.31), (28952.82 28557.68), (…
3 Outram 255 ((28290.63 29757.42), (28596.68 29346.08), (…
4 Rochor 214 ((29383.5 31768.77), (29521.51 31404.98), (2…
5 Novena 182 ((28059.14 33734.79), (28062.48 33523.59), (…
6 Queenstown 176 ((20006.09 30435.35), (20082.89 30453.04), (…
7 Bukit Merah 163 ((24405.42 28144.15), (24581.26 28085.54), (…
8 River Valley 157 ((27291.23 31642.72), (27473.75 31495.65), (…
9 Geylang 145 ((32942.57 32745.19), (33084.45 32747.7), (3…
10 Bedok 120 ((36040.87 33769.19), (36077.61 33347.9), (3…
# … with 33 more rowskallang has the most airbnbs
K cross function
NetSPPA
kallang airbnbs:
kallang_airbnb_sf <- filter(airbnb_sf, neighbourhood == "Kallang")
kallang_airbnb_sfSimple feature collection with 353 features and 4 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 29508.14 ymin: 31119.73 xmax: 33791.74 ymax: 34536.48
Projected CRS: SVY21 / Singapore TM
# A tibble: 353 × 5
name neigh…¹ room_…² price geometry
* <chr> <chr> <chr> <dbl> <POINT [m]>
1 EnsuiteA w Balcony view/acce… Kallang Privat… 57 (31346.67 32390.22)
2 Ensuite M1B within Heritage … Kallang Privat… 50 (30798 32708.67)
3 EnsuiteK w Balcony view/acce… Kallang Privat… 57 (31239.83 32106.04)
4 Heritage apartment: Master r… Kallang Privat… 69 (30505.31 32834.72)
5 Lavender QueenBed *rmFr: 5m … Kallang Privat… 59 (31368.93 32169.07)
6 Life Impact Coaching Kallang Entire… 130 (32498.52 32598.12)
7 1-Pax Small Private Room, Sh… Kallang Privat… 78 (31000.55 32666.66)
8 2-pax Small Private Room, Sh… Kallang Privat… 110 (31000.55 32666.66)
9 1-Pax Small Private Room, Sh… Kallang Privat… 78 (31000.55 32666.66)
10 1950s cozy window singlebed … Kallang Privat… 49 (31207.55 32608.05)
# … with 343 more rows, and abbreviated variable names ¹neighbourhood,
# ²room_typeplot(network_sf)
kallang airbnb plot
tmap_mode('plot')
tm_shape(kallang_airbnb_sf) +
tm_dots() +
tm_shape(network_sf) +
tm_lines()
all airbnbs plot
tmap_mode('view')
tm_shape(airbnb_sf) +
tm_dots() +
tm_shape(network_sf) +
tm_lines()